
library(rio)
library(ggplot2)
library(texreg)
library(tidyverse)
library(AER)
library(gridExtra)

##set the working directory for the replication archive
working_dir <- "~/Dropbox/dueling project/2nd_paper/replication archive/"

source(paste0(working_dir,"Auxiliary Scripts/RegressionRunner.R"))

##load in the main datasets
old_news <- readRDS(paste0(working_dir,"Data/newspapers_historical.RDS"))
all_together <- readRDS(paste0(working_dir,"Data/contemporary_social_capital_data_with_news.RDS"))

old_news$all_papers1884 <- old_news$all_papers

all_together$countyfips <- all_together$fips
all_together <- merge(all_together, old_news %>% dplyr::select(countyfips, all_papers1884), all.x = TRUE)

# all_together$all_papers <- all_together$all_papers1884

##prepare the regressions by defining the names of the dependent variables
varnames <- c("AllOrgs","Putnam","pvote2012","nccs2014",
              "DRUGTOT", "MURDER", "RAPE", "GRNDTOT",
              "Mortality", "Premature mortality rate",
              "% without health insurance"
)

varnames_rev <- c("Num. Social Orgs","Num. Civic Orgs","Voter Turnout (2012)","Num. Nonprofits",
                  "Drug Arrest Rate", "Murder Arrest Rate", "Rape Arrest Rate", "All Arrest Rate",
                  "Mortality Rate", "Premature Mortality Rate",
                  "% without Health Insurance"
)


##run the regressions for both the pretreatment and post treatment cases
# all_est <- as.data.frame(do.call(rbind, lapply(varnames, reg_runner_news, dataset = all_together)))
# all_est$DV <- as.factor(varnames_rev)

##make coefficient plots
# sc_plots_rev <-
#   all_est %>%
#   ggplot(aes(reorder(DV,Estimate), Estimate)) +
#   geom_point(size = 5, position=position_dodge(width=0.5)) +
#   geom_errorbar(
#     aes(ymin = Estimate - 1.65*`Std. Error`, ymax = Estimate + 1.65*`Std. Error`),
#     width = 0.01,
#     linetype = "solid",
#     position=position_dodge(width=0.5)) +
#   xlab("") +
#   ylab("Coefficient on # of Newspapers (logged, 2000)") +
#   geom_hline(yintercept = 0, linetype = "dotted") +
#   coord_flip() + 
#   theme_bw() +
#   theme(legend.position = "bottom", 
#         legend.title = element_text(size = 25),
#         axis.text.x = element_text(size = 25), axis.text.y = element_text(size = 25),
#         legend.text = element_text(size = 25),
#         axis.title.x = element_text(size = 25),
#         axis.title.y = element_text(size = 25),
#         strip.text.x = element_text(size = 25)) +
#   theme(legend.position="bottom")


all_est <- as.data.frame(do.call(rbind, lapply(varnames, reg_runner_news_1884, dataset = all_together)))
all_est$DV <- as.factor(varnames_rev)

##make coefficient plots
sc_plots_rev1 <-
  all_est %>%
  ggplot(aes(reorder(DV,Estimate), Estimate)) +
  geom_point(size = 5, position=position_dodge(width=0.5)) +
  geom_errorbar(
    aes(ymin = Estimate - 1.65*`Std. Error`, ymax = Estimate + 1.65*`Std. Error`),
    width = 0.01,
    linetype = "solid",
    position=position_dodge(width=0.5)) +
  xlab("") +
  ylab("Coefficient on Number of\nNewspapers per 1,000 Residents (1884)") +
  geom_hline(yintercept = 0, linetype = "dotted") +
  coord_flip() + 
  theme_bw() +
  theme(legend.position = "bottom", 
        legend.title = element_text(size = 25),
        axis.text.x = element_text(size = 25), axis.text.y = element_text(size = 25),
        legend.text = element_text(size = 25),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25),
        strip.text.x = element_text(size = 25)) +
  theme(legend.position="bottom")


all_est <- as.data.frame(do.call(rbind, lapply(varnames, reg_runner_news_change, dataset = all_together)))
all_est$DV <- as.factor(varnames_rev)

##make coefficient plots
sc_plots_rev2 <-
  all_est %>%
  ggplot(aes(reorder(DV,Estimate), Estimate)) +
  geom_point(size = 5, position=position_dodge(width=0.5)) +
  geom_errorbar(
    aes(ymin = Estimate - 1.65*`Std. Error`, ymax = Estimate + 1.65*`Std. Error`),
    width = 0.01,
    linetype = "solid",
    position=position_dodge(width=0.5)) +
  xlab("") +
  ylab("Coefficient on the Change in the Number of\nNewspapers per 1,000 Residents (1884-2000)") +
  geom_hline(yintercept = 0, linetype = "dotted") +
  coord_flip() + 
  theme_bw() +
  theme(legend.position = "bottom", 
        legend.title = element_text(size = 25),
        axis.text.x = element_text(size = 25), axis.text.y = element_text(size = 25),
        legend.text = element_text(size = 25),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25),
        strip.text.x = element_text(size = 25)) +
  theme(legend.position="bottom")




ggsave(sc_plots_rev1, 
       file  = paste0(working_dir,"Replication Plots/social_capital_plots_newspapers_1884.pdf"), 
       width = 13, height = 8)


ggsave(sc_plots_rev2, 
       file  = paste0(working_dir,"Replication Plots/social_capital_plots_newspapers_change.pdf"), 
       width = 13, height = 8)


big_plot <- grid.arrange(sc_plots_rev1,sc_plots_rev2)

ggsave(big_plot, 
       file = paste0(working_dir,"Replication Plots/social_capital_plots_newspapers_SI.pdf"), 
       width = 14, height = 10)

